Climate change vulnerability and conservation strategies for tertiary relict tree species: Insights from landscape genomics of Taxus cuspidata

Abstract The unprecedented habitat fragmentation or loss has threatened the existence of many species. Therefore, it is essential to understand whether and how these species can pace with the environmental changes. Recent advantages in landscape genomics enabled us to identify molecular signatures of adaptation and predict how populations will respond to changing environments, providing new insights into the conservation of species. Here, we investigated the pattern of neutral and putative adaptive genetic variation and its response to changing environments in a tertiary relict tree species, Taxus cuspidata Sieb. et Zucc, which is distributed in northeast China and adjacent regions. We investigated the pattern of genetic diversity and differentiation using restriction site‐associated DNA sequencing (RAD‐seq) and seven nuclear microsatellites (nSSRs) datasets. We further explored the endangered mechanism, predicted its vulnerability in the future, and provided guidelines for the conservation and management of this species. RAD‐seq identified 16,087 single nucleotide polymorphisms (SNPs) in natural populations. Both the SNPs and nSSRs datasets showed high levels of genetic diversity and low genetic differentiation in T. cuspidata. Outlier detection by F ST outlier analysis and genotype‐environment associations (GEAs) revealed 598 outlier SNPs as putative adaptive SNPs. Linear redundancy analysis (RDA) and nonlinear gradient forest (GF) showed that the contribution of climate to genetic variation was greater than that of geography, and precipitation played an important role in putative adaptive genetic variation. Furthermore, the genetic offset and risk of non‐adaptedness (RONA) suggested that the species at the northeast edge may be more vulnerable in the future. These results suggest that although the species has maintained high current genetic diversity in the face of recent habitat loss and fragmentation, future climate change is likely to threaten the survival of the species. Temperature (Bio03) and precipitation (Prec05) variables can be potentially used as predictors of response of T. cuspidata under future climate. Together, this study provides a theoretical framework for conservation and management strategies for wildlife species in the context of future climate change.

genotype-environment associations (GEAs) revealed 598 outlier SNPs as putative adaptive SNPs.Linear redundancy analysis (RDA) and nonlinear gradient forest (GF) showed that the contribution of climate to genetic variation was greater than that of geography, and precipitation played an important role in putative adaptive genetic variation.Furthermore, the genetic offset and risk of non-adaptedness (RONA) suggested that the species at the northeast edge may be more vulnerable in the future.These results suggest that although the species has maintained high current genetic diversity in the face of recent habitat loss and fragmentation, future climate change is likely to threaten the survival of the species.Temperature (Bio03) and precipitation (Prec05) variables can be potentially used as predictors of

| INTRODUC TI ON
Habitat loss or fragmentation due to climate change or human activity is recognized as a major threat to global biodiversity (Pykälä, 2019;Rands et al., 2010).Among the various forms of biodiversity, genetic diversity is particularly important in conservation biology, as it provides the raw material for evolutionary change and thus has the potential to adapt to changing environments (Laikre et al., 2010;Nielsen et al., 2022;O'Brien et al., 2022).Normally, genetic diversity might decrease because of increased random genetic drift, inbreeding, and reductions in gene flow due to habitat loss or fragmentation of species (Aguilar et al., 2008;Chung et al., 2023;Lowe et al., 2005;Young et al., 1996).However, habitat loss or fragmentation can also lead to highly specific genetic consequences (e.g., increased genetic divergence, genetic bottleneck and reduced effective population size) for plants because of their different life histories, particularly for tree species with long generation times (e.g., González et al., 2020;Petit & Hampe, 2006;Piotti, 2009).
Forest tree species are largely undomesticated and exhibit local adaptation across heterogeneous environments (Bonan, 2008;Litkowiec et al., 2018).Neutral molecular markers, such as microsatellite genotyping, chloroplast DNA (cpDNA) and mitochondrial DNA (mtDNA) fragments, have been used to reveal the genetic dynamics of forest tree species especially endangered tree species.
In general, low genetic diversity and strong genetic differentiation were revealed in endangered or critically endangered tree species with small population size because of the enhanced genetic drift effect and inbreeding (e.g., Frankham, 2005;Yang et al., 2020), such as T. yunnanensis (Miao et al., 2016), Abies ziyuanensis (Tang et al., 2008), Cathaya argyrophylla (Wang & Ge, 2006).However, high levels of genetic diversity have also been detected in some endangered tree species with current small or scattered population size (e.g., de Sousa et al., 2020;Gitzendanner et al., 2012;Yousefzadeh et al., 2018).
To date, neutral molecular markers are the most common tools for conservation genetic studies, however, genetic variation at neutral loci cannot provide direct information on adaptive selective processes involving the interactions between individuals and their environment (Hoffmann & Willi, 2008;Holderegger et al., 2006).
Recently, with the advent of high-throughput sequencing technologies, adaptive genetic markers have become available even for tree species with large and complex genomes (Feng & Du, 2022;Mackay et al., 2012;Neale & Kremer, 2011).In addition, emerging landscape genomics aims to detect adaptive variation along the landscape by integrating geographic and environmental information to provide unprecedented insights into the evolutionary mechanisms of local adaptation (Borrell et al., 2020;Du et al., 2020;Sork et al., 2013;Wang et al., 2021).F ST outlier analysis and genotype-environment associations (GEAs) in landscape genomics have been suggested as promising tools for detecting adaptive genetic variation in conservation practices (Ahrens et al., 2018;Chung et al., 2023;Schoville et al., 2012).In recent years researchers have increasingly employed these methods to study how spatial heterogeneity promotes population genetic dynamics, and how genomic variation promotes adaptive evolution in forest tree species, particularly conifers (Jia et al., 2020;Mayol et al., 2020;Zhao et al., 2020).
Taxus cuspidata Sieb.et Zucc., a tertiary relict tree species, is a keystone-dominant coniferous long-lived wind-pollinated tree species with a discontinuous distribution in Japan, Korea, Northeast (NE) China and Far Eastern Russia (Kunashir Island, Sakhalin and Primorye) (Chung et al., 1999;Kitamura & Murata, 1987).In China, the species is mainly located at altitudes ranging from 600 m to 1200 m in the Changbai Mountains and neighboring areas.The number of natural populations of T. cuspidata has drastically decreased in the last century, partially because of human disturbance.Therefore, T. cuspidata is considered as "Plant Species with Extremely Small Populations (PSESP)" in China (Wade et al., 2016), however, it was a Least Concern (LC) species in the International Union for Conservation of Nature (IUCN) (Katsuki & Luscombe, 2013).Previous studies based on paternal inherited mtDNA and cpDNA fragments have shown high degrees of genetic diversity in T. cuspidata (Cheng et al., 2015;Kozyrenko et al., 2017;Su et al., 2018).However, the contributions of geography and climate to genetic variation in T. cuspidata remain unclear.
Here, we sampled T. cuspidata natural populations throughout the species' distribution range and combined their neutral and nonneutral genetic variation using nSSRs and SNPs.We aimed to (1) es-   S1.Leaf samples were labeled and stored in plastic bags with silica gel for DNA extraction.Voucher specimens of each population were deposited at Beijing Forestry University.

| DNA isolation and microsatellite genotyping
Total genomic DNA was extracted from leaf tissue using Plant Genomic DNA Kit (Tiangen, Beijing, China) following the manufacturer's protocol.The DNA quality was initially checked by 1% agarose electrophoresis gel and then the concentration was measured by an ultramicro-spectrophotometer (Thermo Fisher, USA).Seven polymorphic nuclear microsatellite (nSSR) loci were used for genotyping all 200 individuals after initial screening of 24 nSSR primers that had been already applicable to T. cuspidata (Cheng et al., 2015;Dubreuil et al., 2008;Kondo, 2016) (see Table S2 for primer details).
The PCR conditions followed Du et al., 2017 and the PCR products were analyzed using an ABI PRISM 3730 Genetic Analyzer (Applied Biosystems, USA).Subsequently, the alleles were scored using GENEMARKER v. 2.2 (Softgenetics, USA) and the genotype was checked twice.

| Double-digest RADseq-derived SNP dataset
A subset of samples (130 out of 200 individuals in 19 populations) was used for SNP dataset collection by a double-digestion restriction fragment-based procedure (ddRAD-seq) (Peterson et al., 2012) (Table S3).The subset of sampling by RAD-seq covered almost the entire distribution range to capture the vast majority of diversity.
For each sample, the genomic DNA was processed for library construction and sequencing.Briefly, DNA was double-digested with restriction enzyme TaqI and MseI, followed by the ligation of Illumina adaptors.Ligation products were size-selected about 500 bp and amplified by kapa HotStart ReadyMix (cat.no.KK2601; KAPA Biosystems) with 13 cycles.Paired-end sequencing (2 × 150 bp) was performed using an Illumina HiSeq2000 platform at Majorbio Pharm Technology Co., Ltd., Shanghai, China.
The data quality control was assessed by FastQC v0.11.7 (Andrews, 2010).We removed adapter sequences and low-quality bases (Phred quality <20) from raw data using Trimmomatic 0.36 (Bolger et al., 2014).Then the reads were all trimmed to 120 bp and designed for SNP genotyping using Stacks v1.48 (Catchen et al., 2011(Catchen et al., , 2013)).In detail, ustacks was used to assemble our sequences into de novo, cstacks was used to build a catalog of loci, sstcaks was used to match loci against the catalog and the Populations was used to output SNPs.A strict criteria was used to filter SNPs using VCFtools (Danecek et al., 2011): first, we only kept individuals that represented at least 60% of the SNPs; second, a filtering parameter of 0.8 was used to avoid the influence of missing data, i.e., more than 80% of individuals in a population were required to process a locus; third, a minor allele frequency (MAF) < 0.1 was used to filter the data and to reduce the likelihood of false-positive results due to spurious correlations.

| Genetic diversity and structure
We estimated genetic diversity including mean effective number of alleles (N E ), mean observed heterozygosity (H O ), mean expected heterozygosity (H E ) and mean number of different alleles (N A ) for each population by GenAlEx 6.5 for nSSR dataset (Peakall & Smouse, 2006).We defined two types for ddRAD-seq dataset: neutral genetic variation based on neutral SNPs and non-neutral (putatively adaptive) genetic variation based on outlier SNPs (i.e., outliers were detected by F ST and GEA analysis; see section outlier detection below).We estimated the mean frequency of the most frequent allele at each locus (P), mean nucleotide diversity (π), H O and H E for all, neutral, F ST -outlier and GEA-outlier SNPs by Populations module in Stacks (Catchen et al., 2013).
We conducted a principal component analysis (PCA) to produce a lower dimensional subspace that captured most of the variation in each of three data types.For nSSRs, "adegenet" (Jombart & Ahmed, 2011) was used for the PCA in R 3.6.1 (R Core Team, 2019).PLINK v.1.07was used for the PCA in all, neutral and outlier SNPs (Purcell et al., 2007;Zheng et al., 2012).
Population structure was also accessed using Bayesian clustering.
For nSSRs, we identified population genetic structure through the Bayesian Markov chain Monte Carlo (MCMC) clustering method implemented in STRUCTURE v 2.3.4 (Pritchard et al., 2000).Twenty independent runs were performed for each value of K (1-10) using 200,000 generations for the MCMC cycles and 100,000 generations for the burn-in by STRUCTURE HARVESTER (Earl & vonHoldt, 2012).The most likely number of clusters (K) was determined using ΔK and LnP(K) statistics, according to Evanno et al. (2005).CLUMPP (Jakobsson & Rosenberg, 2007) and DISTRUCT v1.1 (Rosenberg, 2004) were used to aligned independent runs and visualize the bar plots of the individual's probabilities of population membership.ADMIXTURE v1.3.0 (Alexander et al., 2009;Alexander & Lange, 2011), which can deal with large data with fast speed, was applied to assess population structure of genetic variation for all, neutral and outlier SNPs.We ran ADMIXTURE with default 5-fold cross-validation (−cv = 5) to select the optimal K from one to ten, the best K exhibits low cross-validation error (CV error) opposed to others (Alexander & Lange, 2011).
We conducted a hierarchical analysis of molecular variance (AMOVA) to quantify the degree of genetic divergence among and within populations for our nSSRs, all SNPs, neutral SNPs, F SToutlier SNPs and GEA-outlier SNPs by Arlequin 3.5 (Excoffier & Lischer, 2010).The significance of genetic differentiation differences was evaluated using 10,000 permutations.

| Climatic variables
Bioclimatic variables of current climate (representative of  and future climate (2070: average for 2061-2080) were downloaded from the WorldClim Version2 (http:// world clim.org/ version2, Fick & Hijmans, 2017) at spatial resolutions of 30 s (~1 km 2 ).A total of 31 climate variables were extracted by longitudes and latitudes of population sites based on raster layers by ArcMap10.2.To avoid biased estimates of model coefficients and spurious significance levels resulting from multicollinearity, we excluded highly correlated climate variables with the threshold values of 0.7 using a variance inflation factor (VIF) test in "usdm" R package (Marquaridt, 1970;Naimi et al., 2014;R Core Team, 2019).Six climate variables were finally retained: isothermality (bio03), maximum temperature of warmest month (bio05), mean temperature of driest quarter (bio09), precipitation seasonality (bio15), precipitation in May (prec05) and precipitation in October (prec10).We also carried out Shapiro-Wilk and Levene test (package car, Fox & Weisberg, 2019) to explore each variable's data distribution and homogeneity.We found that the data for each variable was not normally distributed nor homogeneous; therefore, we performed the Kruskal-Wallis rank sum test to explore the significance level for each variable.

| Outlier detection
Three algorithms were used to detect potentially adaptive loci, including one F ST outlier analysis approach and two genotypeenvironment associations (GEAs) outlier detection approaches.
For the F ST outlier analysis, we use Bayesian approach in BayeScan 2.0 to directly estimate the posterior probability of a given locus under selection (Foll & Gaggiotti, 2008).We used the following parameter values: sample size of 5000, 20 pilot runs with 5000 run length, 50,000 burn-in iterations and thinning interval of 10. Prior odds for neutral model were set to 10 and SNPs with q < 0.01 (− log 10 q > 2) were considered as F ST -outliers (Puebla et al., 2014).
For genotype-environment associations (GEAs) outlier detection approaches, we firstly applied latent factor mixed models (LFMM) by package LEA in R (R Core Team, 2019).We use a hierarchical Bayesian mixed modelling approach to identify allele-environment correlations, while modelling residual population structure with 'latent factors' (Frichot et al., 2013;Frichot & François, 2015).In LFMM, environmental variables were tested separately and introduced into each model as fixed effects, and the number of latent factors (K = 3) was selected to account for genetic structure by sparse nonnegative matrix factorization (SNMF) (Figure S7).The parameters were set as follows: 100,000 iterations with a burn-in of 50,000 iterations and ten replicate runs.Significant outliers were determined as SNPs with p-values of p < 10 −3 or -log 10 (p-value) > 3. We secondly used Bayesian generalized linear mixed models (BayEnv) to detect the correlation between SNPs and environmental variables with a neutral dataset generated by BayeScan as a null model (Günther & Coop, 2013).We initially computed a null covariance matrix of relatedness between populations, over 100,000 iterations and five independent runs.We then tested all SNPs (including those initially identified by BayeScan) under an alternative model where allele frequencies are determined by a combination of the covariance matrix and an environmental variable.We performed our analysis independently across six environmental variables, and using Bayes factor (BF) to evaluate the posterior probability that each SNP is under selection across independent environmental variable.We also performed nonparametric Spearman's rank correlation as alternative tests to the BF and detect the correlation between ranks of SNP allele frequencies and environmental variables.Therefore, we considered the SNPs in the top 1% of BF values (BF > 3) and top 10% of the absolute value of Spearman rank correlation coefficients (ρ) as significant outliers.

| Multivariate relationship between genetic variation and environmental gradients
To estimate how the genomic variation is influenced by climate or geographic variables, we firstly performed redundancy analyses (RDAs) and partial redundancy analyses (pRDAs) to detect linear relationships between genetic variations and multivariate climatic gradients, using "vegan" R package (Oksanen et al., 2019).We constructed two pRDAs models to differentiate the independent effects of climate and geography by reciprocally constraining one of the two factors.Briefly, climatic effects were conditioned on the effects of geography (pRDAenv to extract pure effects of the climate) and vice versa (pRDAgeo for pure effects of geography).All SNPs, neutral and non-neutral genetic variation with the six climate variables and geographic variables (longitude and latitude) were considered as predictor variables.Statistical significance was evaluated from 999 permutations.
We then employed gradient forest (GF) analyses to explore the nonlinear relationships between environmental variables and their contribution to genetic variation using R package GradientForest (Ellis et al., 2012).Gradient forest is a machine learning, regression tree approach that allows for exploration of nonlinear associations of spatial, environmental and allelic variables (Ellis et al., 2012;Fitzpatrick & Keller, 2015).We conducted a GF analysis on six climate variables to assess the relative importance of each predictor variable using weighted R 2 value, split importance; Ellis et al. (2012).Split importance, a measure of the amount of variation explained, is high in positions along the gradient where allelic change is large.We ran a gradient forest with 500 regression trees per SNP, maxLevel = log 2 (0.368n)/2, and variable correlation threshold of 0.5 to calculate conditional variable importance as recommended (Ellis et al., 2012;Fitzpatrick & Keller, 2015).The default values of the other parameters were carried out for each GF model.

| Adaptation potential of Taxus cuspidata
We employed two different methods to estimate the adaptive potential of T. cuspidata.We selected the Global Climate Model BCC-CSM1-1 (IPCC, 2014) under two contrasting representative concentration pathways (RCPs), including a low-emission scenario (RCP 2.6) and a high-emission scenario (RCP 8.5) in 2070 years for future climate scenarios.We first used GF to estimate the genetic offset based on all SNPs, F ST -outlier SNPs and GEA-outlier SNPs under future climatic conditions.Genetic offset, a measure of the magnitude of genetic change required between present and future climate (Fitzpatrick & Keller, 2015) was used to predict genetic variation across grid cells of NE China by "GradientForest" package (Ellis et al., 2012) in R (R Core Team, 2019).For each grid cell, the Euclidian distances between the current and future genetic compositions were calculated and served as the metric for genetic offset (Ellis et al., 2012).
We also performed a risk of non-adaptedness (RONA) to pre-

| RAD-based SNP data
We attained a total of 373,894,701 loci with the mean depth of coverage for filtered SNPs at 4.89 (range: 3.58-7.47)and mean proportion of loci with missing data per sample is 0.46 (range: 0.07-0.98)(Table S3).A total of 37 individuals with high proportion missing data were discarded, and the retaining 93 individuals is at least four individuals per population (Table 1).Stacks initially recovered 144,575 SNPs and 16,087 high-quality SNPs were retained after stringent quality control.

| Genetic diversity and differentiation based on nSSRs and SNP datasets
For nSSRs, observed and expected heterozygosity estimated per population ranged from 0.26 to 0.46 and 0.28 to 0.48, respectively (Table 1, Table S5).The expected heterozygosity (H E ) and nucleotide diversity (π) ranged from 0.26 to 0.30 and 0.29 to 0.33 in all (Table S4) and neutral SNPs (Table 1).For F ST -outlier SNPs, H E and π ranged from 0.17 to 0.38 and 0.18 to 0.40, and for GEA-outlier SNPs, H E and π ranged from 0.28 to 0.34 and 0.31 to 0.36 (Table 1).S1, S2) based on both nSSRs and RAD-seq datasets were consistent with Bayesian clustering analysis, which failed to detect clear population genetic structure in all populations.The most likely number of clusters (K) was two for the nSSRs (Figure 1a, Figures S3, S4), one for all and neutral SNPs (Figure 5b, S5a), and three for outlier SNPs (Figure 1b, Figures S5c, S6).These results showed a mixed genetic makeup of clusters in all populations.

PCA (Figures
The results of AMOVA indicated that a large proportion of genetic variation occurred within populations, with low levels of genetic differentiation (F ST = 0.07) in nSSRs, all SNPs and neutral SNPs; For outlier SNPs, high F ST values (F ST = 0.32) in F ST -outlier SNPs and moderate F ST values (F ST = 0.1) in GEA-outlier SNPs (Table 2).

| Outlier detection and environmental association
We identified 63 SNPs as putative outliers and did not detect any significantly low outlier F ST values that would be indicative of balancing or purifying selection using BayeScan (Figure S8).According to the significance test, each climate variable was significantly different among the populations (Kruskal-Wallis chi-squared = 92, p < 0.001).We identified 279 and 286 putative adaptive SNPs that were significantly associated with at least one climatic variable using LFMM (Figure S9) and BayEnv (Figure S10), respectively, with five common SNPs (Figure 2c).More putative adaptive SNPs were significantly associated with precipitation variables than temperature variables (Table S6).A total of 598 outlier SNPs were identified as non-neutral genetic variations and 15,489 neutral SNPs as neutral genetic variations.

| Environmental associations with genetic variation
RDA revealed that a large proportion of genetic variation among populations was associated with the six climate variables.The RDA results were similar to that of pRDA (Tables S7 and S8; Figures S11a,   S12); here, we explained the results of pRDA (Table 3, Figures S11b,   S13).For nSSRs, genetic variation was found to be mainly associated with climate variables: 5% of the nSSRs variance was explained by climate, while geography only accounted for 1% (Table 3).For all SNPs and neutral SNPs, 10% of the explained variance was explained by climate and 3% by geography (Table 3, Table S8).For outlier SNPs, TA B L E 1 Genetic diversity of Taxus cuspidata based on neutral and non-neutral genetic variation.21% of the F ST -outlier SNPs variance was explained by climate and 4% by geography; 13% of the GEA-outlier SNPs variance was explained by climate and 3% by geography (Table 3).

Non
For all SNPs, the GF indicated that isothermality (bio03) was the most important climate response factor (Figure 2a).However, for outlier SNPs, precipitation in May (prec05) and precipitation seasonality (bio15) were the most important climate response factors, respectively (Figure 2b,c).

| Prediction of population vulnerability under future climate change
In the three GF models by all SNPs, F ST -outlier and GEA-outlier SNPs, the predicted turnover in allele frequencies across the landscape followed a southwest to northeast direction: northeastern populations are expected to have higher genetic offsets in the Changbai Mountains and adjacent areas in the future (Figure 3a-f), the degree  S1.

Source of variation df
Percentage of variation (%)  of genetic offset is slightly higher under RCP 8.5 than that of RCP 2.6 (Table S9).The RONA revealed that the most represented environmental variables were prec05, bio03, and prec10 (Table S10).
Under the combined effect of the three most represented climate variables, population THL showed the highest RONA value for RCP 2.6 (Figure 4).MDJS, HCM and YBJ populations showed the highest RONA value for RCP 8.5 in 2070 (Figure 4).The THL and MDJS populations had a lower adaptive potential for prec10 by 2070.

| DISCUSS ION
Our study integrated population genetics and landscape genomic methods to explore neutral and non-neutral genetic variations in T. cuspidata populations in the Changbai Mountains.The neutral genetic variation dataset revealed high genetic diversity and low genetic differentiation, whereas non-neutral genetic variation showed high genetic differentiation.GEAs revealed that the contribution of climate to genetic variation was greater than that of geography, and precipitation played an important role in the putative adaptive genetic variation in response to climate change.We found that populations in the northeast range will be more vulnerable to future climate change as suggested by putative adaptive genetic variation.
Our study provides insight into how neutral and putative adaptive genetic variation interacts with the environment, which is essential for future conservation and management of natural populations.

| High genetic diversity was maintained in present fragmented populations
We integrated the patterns of neutral and non-neutral genetic variation in T. cuspidata using nSSRs and RAD-seq datasets.In this TA B L E 3 Summary and partitioning of the variance associated with climate and geographical variables based on pRDA in neutral and nonneutral genetic variation.S9.  et al., 2017), cpDNA and mtDNA fragments in China (Su et al., 2018), and were similar to that of T. baccata from Poland by microsatellites (Litkowiec et al., 2018).In addition, admixture analysis revealed substantial gene flow among the populations.For fragmented and isolated T. cuspidata populations, substantial gene flow may compensate for the genetic barrier caused by fragmented habitats (Browne & Karubian, 2018;Hu et al., 2008;Petit et al., 2005).
Theoretically, fragmentation of species usually has a negative effect on genetic diversity compared to intact populations (Nybom, 2004;Vranckx et al., 2012).One possibility for the high genetic diversity and gene flow might be the outcrossing nature of the species.T. cuspidata is a wind pollination and typical birddispersed plant.Gene flow may be accomplished by both seeds and pollen, and pollen often disperses over longer distances than seeds (Jordano, 2017;Kremer et al., 2012).However, this seems to be difficult, particularly for long-distance intervals among the studied populations.For yew, a multi-model approach indicated that 95% of the seed and pollen dispersal distances were <109 m and 704 m, respectively (Chybicki & Oleksa, 2018).An alternative explanation for high genetic diversity is the delayed sexual maturity of the species.The unbalanced ratio of females to males of T. cuspidata (1:2.3) was found in this region (Long et al., 2021).The sex ratio in dioecious plants might significantly affect genetic diversity (Rosche et al., 2018), e.g., the dioecy may increase genetic diversity due to obligate outcrossing, which has been found in Pherospheara hookeriana (Worth et al., 2021).Isolation occurred mainly throughout the last few decades and was probably not long enough to impact current genetic diversity (Münzbergová et al., 2013;Rosche et al., 2018).Therefore, the short duration of fragmentation tends to maintain the standing variations in long-lived plants, which may delay adverse effect in the progenies of T. cuspidata (Aguilar et al., 2008;Vranckx et al., 2012).

| Signature of natural selection in T. cuspidata
Outlier detection test is an effective approach to identify the presence of putatively adaptive genetic variation.Although outlier detection tests can produce false outliers due to confounding factors (Schoville et al., 2012), tests with various demographic hypotheses can be utilized and compared to solve this issue, and common loci detected in consensus are more likely to be actual targets of selection (Ahrens et al., 2018).In this study, we applied three outlier methods to detect different sets of outlier SNPs that showed signatures of selection.Temperature is a key factor influencing the growth and phenology of tree species (Schmiege et al., 2021;Vitasse et al., 2009).RDA and GF indicated that temperature was found to be associated with our nSSRs, neutral SNPs and all SNPs sets.An explanation to this temperature-driven genetic variation is that temperature could affect the connection between populations (i.e., gene flow), which has been observed in yews in the northern Italy (Mercuri et al., 2013).Yews are also highly sensitive to extreme temperature changes, including T. globosa (Antúnez, 2021) and T.
baccata (Mayol et al., 2020;Moir, 1999).Low temperatures can encourage the rooting of yews during the early growth stage (Muñoz-Gutiérrez et al., 2009).Therefore, the maximum temperature of warmest month (bio05) and the mean temperature of driest quarter (bio09) might affect the growth and distribution of T. cuspidata.
We also found that precipitation was an important variable associated with putative adaptive genetic variation.Water availability may be an important selection factor for yew species, especially during the last interglacial period, as indicated by the association study in the common yew, T. baccata (Mayol et al., 2020).A previous study using ecological niche modeling has also found that the best suitable distribution time for T. cuspidata is during the last glacial maximum and the most restricted distribution time for the species is during the last interglacial period (c.130 ka), suggesting the coldtolerant and wet-sensitive features of the species (Su et al., 2018).
Therefore, water availability will be a challenge for T. cuspidata when facing the changing environments, as suggested by other yew species (Linares, 2013).

| Species vulnerability and conservation strategies
Assessing the vulnerability of species to climate change can provide insight into the potential risk of species persistence in future climate scenarios (Capblancq et al., 2020;Feng & Du, 2022).
Genetic prediction of adaptation to future climate indicated that populations (e.g., MDJHP, YJX and MDJS) in the northeast area (Heilongjiang and Jilin Provinces) exhibited a higher risk of genetic maladaptation than other populations.These populations showed low potential to adapt to changing environments in small or isolated habitats.
For GF genetic offset, the predictive power of genomic offset estimates on fitness effects is increasingly being assessed through experimental and simulation studies, showing promising results (Láruson et al., 2022).The GF offset results in our study differed not much between the two future climate scenarios.An explanation for the phenomenon is that T. cuspidata has a long lifecycle with long generation intervals, in which the rates of emergence and spread of novel adaptive alleles in populations through de novo mutations are likely to be too slow to respond to climate change.
For all SNPs and putatively adaptive SNPs, northeastern populations showed higher genetic offsets and vulnerabilities than the other populations.Su et al. (2018) found that ecological niche modeling showed a contraction trend in the distribution of T. cuspidata by 2070, with northeastern populations located in a contracting area.Láruson et al. (2022) showed that higher genetic offset may be caused by genetics drift at small population rather than a selectiondriven response.Negative associations between GF offset and population size have also been found in previous bird studies (Bay et al., 2018;Ruegg et al., 2018).Therefore, the high genetic offset of T. cuspidata may also be caused by natural selection and/or small population size.
Based on our study and previous reports, we came to the conclusion that T. cuspidata is at risk in China with number of challenges, such as fragmentation and habitat loss, climate change, disturbance and the unbalanced ratio of females to males (Long et al., 2021;Su et al., 2018;Wang et al., 2019).Our study specifically recommends adopting in situ conservation strategies for THL, LJD and YBHSP populations with high adaptive genetic diversity but vulnerable in the future.In addition, populations displaying low adaptive variation, including MDJHCH, MDJHP and YBJ, could benefit from ex situ conservation strategies.Furthermore, for small or isolated populations characterized by limited genetic diversity and high vulnerability (e.g., BSSCZ, YBJ), genetic rescue and assisted gene flow techniques should be considered to facilitate their adaptation to climate change.
It is worth noting that our study represents a pioneering effort to integrate neutral and adaptive genetic variation to the conservation of a threatened tree species.However, to gain deeper insights into the molecular mechanisms underlying the threats faced by T.
cuspidata, future research endeavors should include comprehensive investigations encompassing whole-genome or transcriptome sequencing, along with the inclusion of phenotypic and fine-scale environmental data.

ACK N OWLED G EM ENTS
timate the genetic diversity and genetic differentiation within and among populations, (2) understand the contribution of geographic and climatic factors to genetic variation and (3) predict future adaptability of the species using landscape genomic tools under a scenario of future climate change.Our study lays the groundwork for understanding the molecular mechanisms of local adaptation, and provides a theoretical foundation for the conservation and management of important conifer tree species in East Asia.
response of T. cuspidata under future climate.Together, this study provides a theoretical framework for conservation and management strategies for wildlife species in the context of future climate change.Field sampling In our study, we collected leaf material in natural populations of T. cuspidata in the Changbai Mountains and adjacent areas.The sampling points were strategically distributed across the area, spanning from the southwest to the northeast of Changbai Mountains and adjacent areas.A total of 200 individuals of T. cuspidata were sampled from 19 natural populations and each population was separated by at least 30 km.The detail sampling information of each population was listed in Table the adaptive potential of T. cuspidata to future local climate using default settings in PYRONA v0.3.6 (Pina-Martins et al., 2019;Rellstab et al., 2016).PYRONA ranked the environment factors by the number of associations and provided an average RONA value weighted by the regression R 2 value for each predictor factors.The "RONA value" here indicates the mean difference between current and future expected allele frequencies, which was used to estimate the expected allele frequencies in future based on the present allele frequencies and environmental variables.The inferred linear model was used to predict expected allele frequencies in future environmental gradients in 2070 according to the Global Climate Model BCC-CSM1-1 (IPCC, 2014) under RCP 2.6 and RCP 8.5.

FI G U R E 1
Geographic distribution and sampling sites of Taxus cuspidata.The pie chart shows the genetic clustering of each population based on (a) neutral genetic variation by microsatellites (SSRs) and (b) non-neutral genetic variation by outlier SNPs.The red line indicates the distribution range of T. cuspidata.Orange, green, blue and yellow represent genetic cluster assignments; the details of the cluster assignments were showed in Table

Note:
Abbreviation: df, degree of freedom.

FI G U R E 3
Genetic offset to future climate change predicted for all (a, d), F ST -outlier (b, e) and GEA-outlier SNPs (c, f) of Taxus cuspidata in 2070.(a-c), RCP 2.6; (d-f), RCP 8.5.Red and green indicate high and low genetic offset, respectively.The details of the genetic offset values were showed in Table study, high current genetic diversity and low genetic differentiation were observed across the extant geographical range in China.The present results are in agreement with previous T. cuspidata genetic studies by cpDNA fragments in Russian and South Korea (Kozyrenko High genetic differentiation was also observed in the putative adaptive genetic variation, the high F ST values may be due to the fixation of different alleles in the local populations.This might also be caused by natural selection in response to changing environments, and adaptive processes may contribute to high F ST values.Therefore, natural populations of T. cuspidata primarily stem from intense selective pressure imposed by environmental factors.

F
I G U R E 4 RONA of Taxus cuspidata under RCP 2.6 (a) and RCP 8.5 (b) prediction models for 2070.Bars represent weighted means (by R 2 value), and grey error bars represent the standard error.4.3 | Genetic variation associated with adaptation to climate change Precipitation and temperature are main climatic factors influencing plant's distribution and survival (Cuervo-Alarcon et al., 2018).

-neutral genetic variation nSSRs Neutral SNPs F ST -outlier SNPs GEA-outlier SNPs
Abbreviations: N, number of individuals used for nSSRs and RAD-seq in parentheses; N A , No. of different alleles and N E , effective number of alleles; H o and H E , observed and expected heterozygosity; P, mean frequency of the most frequent allele at each locus; π, mean nucleotide diversity.| 7 of 15 LUO et al.